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Abstract. We investigate the ground state properties of a polarized dipolar Bose- 
Einstein condensate trapped in a triple-well potential. By solving the dipolar Gross- 
Pitaevskii equation numerically for different geometries we identify states which reveal 
the non-local character of the interaction. Depending on the strength of the contact 
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. and dipolar interaction we depict the stable and unstable regions in parameter space. 
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1. Introduction 

The physics of cold atoms in optical lattices is an active field of research in both 
experiment and theory [H [2j [3j HI [5]. Interaction between the atoms and tunneling 
across the lattice are easily controlled using Feshbach resonances and by tuning the 
intensity of the external lasers, respectively. Already in one of the first experimental 
realizations of the system a quantum phase transition between a Mott insulator and a 
superfluid has been shown [3J. Together with the ultra-precise spatial resolution [6j [7] 
this system is a good candidate for quantum simulators j8] or, in the far future, even an 
element of a new generation of computers [9]. 

The search for new phases is ongoing and dipolar interactions present additional 
possibilities [TOl [TT| [12]. The major feature of the dipolar interaction is its long range 
character. Thus, the new phases are expected to reveal inter-site effects even in the 
case of suppressed tunneling. The first indication of such a phenomenon has already 
been shown in the dynamical properties of a Bose-Einstein condensation of very weakly 
interacting 39 K [13] and in the study of the stability of 52 Cr, loaded in a ID optical 
lattice [H] . Furthermore, the anisotropy of the dipolar interaction allows to control both 
inter- and on-site interactions by choosing the appropriate geometry of the lattice sites 
with respect to the polarization direction of the dipoles [15]. The easiest model systems 
consist of a few linked wells. In particular, the double-well system has received a lot of 
attention [161 [T7J. As the entanglement between the two macroscopically occupied modes 
has been demonstrated, it may be considered as an extension of a qubit [HJ [19]. On 
the other hand, many interesting effects, including Josephson oscillations and quantum 
self trapping, were observed in the frame of the mean field approximation |20j . 

In a recent discussion, the triple-well potential, loaded with a dipolar gas, was 
studied |21j. An extended Bose- Hubbard model is used to describe the system, as in 
most related references concerning optical lattices [TQl [23] . This model assumes fixed, 
occupation-independent parameters. For increasing particle numbers and interaction 
strength, this approximation is less reliable. Due to the interaction, the on-site spatial 
distribution of the atoms shrinks in the case of an attractive gas and broadens if the 
interactions are repulsive. The most dramatic case occurs for a dipolar gas when the 
number of atoms or the strength of interaction is above a critical value. The sample 
collapses and then explodes in a so called Bose-Nova [23]. As the shape of the atomic 
cloud changes with the interaction, the parameters of the Bose-Hubbard model cannot 
be uniquely defined. In this paper we study the triple-well case using a mean field 
approach. Within this picture, the ground states for different geometries and interaction 
strength are discussed. We consider experimentally relevant parameters, especially for 
magnetic dipolar gases like Cr or Dy. 

The paper is organized as follows. In section [2] we describe both the Bose-Hubbard 
model and the mean field approach for a dipolar gas in an external triple-well potential, 
which is modelled by overlapping Gaussian wells. We discuss the limitations of the Bose- 
Hubbard approach and then switch to the mean field picture. An important new aspect 
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Figure 1. (a) Cut of the triple- well potential along symmetry axis. It is modelled by 
three overlapping Gaussians of width w. (b) Repulsive and (c) attractive configuration 
of the dipoles. 



arises. Depending on the geometry and the interaction parameters, the ground state 
solution may be unstable. In section [3] we present a phase diagram for two specifically 
chosen geometries. The mean field results are compared to the ground states computed 
with the Bose-Hubbard model. The interesting states revealing the role of the inter-site 
effects are identified. 



2. System and model 

We consider a Bose-Einstein condensate consisting of iV dipolar particles with a dipole 
moment d which may be of either electric or magnetic origin. A strong external field 
is orienting the dipoles such that they all point in the same direction. The particles 
are subject to an external potential V(r) with three (nearly) equivalent minima, see 



figure 1(a), which is modelled by three overlapping Gaussian wells: 
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The parameter Vo determines the centre depth of the individual wells and the widths Wi 
parametrise the geometry of a single well (size in each direction). The spacing between 
the wells is given by I. A potential like this can be created by means of focused Gaussian 
laser beams [2T] . 

The particles are interacting via contact and dipolar interactions. The contact 
interaction is fully characterized by the scattering length a. The dipolar interaction 
between polarized dipoles at positions r 1 , r 2 is given by 

1 -3cos 2 ($) 
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where r = \r\ = \r\ — r 2 \ is the inter-particle distance and i? is the angle between r and 
the dipole moment d. The factor dd is equal to d 2 Ho/47r for magnetic, and d 2 /47reo 
for electric dipoles. In analogy to the scattering length one introduces the length scale 
Odd = m Cdd/3h 2 , characterizing the strength of the dipolar interaction [25J. 

Throughout this work we are using a dimensionless system by measuring all lengths 
in units of the spacing I, all energies in units of h 2 /ml 2 and time in units of ml 2 fK (m 
is the mass of the dipolar particles). We will keep the same notation for quantities with 
and without units, though. 



2.1. Extended Bose-Hubbard model 



Following the discussion in [21J, we outline how a Bose-Hubbard model can be derived for 
the dipolar triple-well system. First we assume that the three minima of the potential 
are well separated, such that the on-site wave functions fa(r) for each site may be 
described by a single function: fa = 0(r — T*j), with V{ the centre of the i-th well. The 
field operator $(r) = $^i=i ^ L i4 > i{ r ) can then be written in terms of the annihilation 
operators a i at site i. Interpreting these operators as annihilating a particle at site i 
is correct as long as the potential is deep enough such that the overlap between the 
wave functions of two adjacent sites is small. The Hamiltonian is then expressed in 
Bose-Hubbard form as 

H = - J a\a 3 + — J^n* (n* - 1) + C/i(n x n 2 + h 2 h 3 ) + U^h^ (3) 

(id) i=l 

where the number operators are defined as hi = a\a i and sums over neighbouring 
sites. The hopping rate is given by 

J = - /d 3 r0*wf-iv 2 +y(r)')0 m (r) (4) 



2 

and the interaction is parametrised by the three parameters 

U k = Jd 3 r J dV |<Mr)| 2 \fa +k (r')\ 2 (^a5(r - r') + V dd (r - r')) . (5) 

The on-site interaction Uq includes parts of both dipolar and contact-interacting origin 
whereas the inter-site couplings U\,U2 only depend on the dipolar interaction, as the 
density-density overlap is negligible. For point-like, tightly localized wave functions fa, 
the inter-site couplings satisfy U 2 = Ui/2 3 as the dipolar interaction falls off like r -3 . 
The following results, however, are also valid for extended wave functions and U 2 = U\/a 
with 4 < a < 8. 

In the special case of J = the model can be solved analytically and four distinct 
phases appear [21J. We quickly review them here to compare with our results. For 
Uq > and U\j |C/o| < 8/15 as well as for U§ < and U\j \Uq\ < —8, the phase A is 
present with 

S(U - U x ) 



ni = n 3 



24U - 31C/i 



-N 



(6) 
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where |_-J denotes the integer part. The other three phases are described by a single fixed 
ground state. Phase B appears for U > and 8/15 < Ui/ \ U \ < 8 and is characterized 
by the ground state n\ = n 3 = N/2. For Uq > and XJ\/ \Uq\ > 8, as well as for Uq < 
and Ui/ \U \ > —1, phase C is present where all particles are occupying a single well. 
The central well is favoured if (weak) tunneling is present and thus we describe this phase 
by n 2 = N. The remaining part of parameter space U < 0, — 8 < U\/ \Uq\ < —1 is filled 
with phase D, having two degenerate states with n x = n 2 = N/2 or n 2 = n 3 = N/2. 

2.2. Restrictions of the Bose-Hubbard model 

There are two main issues with the Bose-Hubbard approach that we will address in this 
section. Both restrictions arise from the assumption that a ground state wave function 
(j) exists which does not depend on the number of particles and the interaction strength. 

The Bose-Hubbard method intrinsically leads to a stable ground state solution 
as this is a premise of the model. This assumption, however, does not hold true in 
general for interacting quantum gases. As we would like to describe both repulsive and 
attractive interactions, this problem is of relevance in our system. The stability issue 
will be discussed in section 12.41 

The second assumption is that the parameters J and Uk, which are calculated by 
means of the single particle wave function 0, are constant for all particle numbers N. 
This approximation is certainly good for small particle numbers and small values of 
a, add- We demonstrate, however, that it is not well suited in our case. 

For simplicity we consider a purely contact interacting Bose-Einstein condensate 
of N particles in a spherically symmetric harmonic well with frequency u>. We 
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Figure 2. Interaction energy of a contact interacting gas in a harmonic trap, compared 
to the energy term Uq/2 ■ N(N — 1) of the Bosc- Hubbard model. The quadratic scaling 
is only reasonable for a/a^N <C 1. The strongly interacting regime is well described 
by the Thomas-Fermi approximation. 
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calculate the interaction energy as a function of a/a^N, where a^ = ^Jn~J 



mu is the 



harmonic oscillator length. The Bose-Hubbard approach suggests a quadratic scaling 
Ei n t = Uq/2N(N — 1) with the number of particles. 

For a/a^oN <C 1 this relation is a valid approximation, see figure [2j However, we 
have in mind a system of at least 100 atoms with a typical scattering length of a ~ 5 nm 
and traps with a width of a^o ~ 500 nm. The resulting factor oia/a^N « 1 is just in the 
crossover region of the diagram. For this value the interaction energy calculated by the 
quadratic term is already 35% off, compared to the numerical simulation. The extension 
to dipolar interacting gases is further increasing the problem. Inter-site repulsion or 
attraction can lead to changes of the neighbouring on-site wave functions. 

2.3. Mean- field approach 

For reasons being apparent now, we will use an alternative approach and describe the 
system in a mean-field picture. We stress that this approach requires, contrarily to the 
Bose-Hubbard model, that the particle number iV 3> 1. The Gross-Pitaevskii equation 
for our system is given by [26J 



where we have introduced the condensate wave function ty(r,t) which we normalize to 
unity. The dipolar interactions are included by the mean-field potential 



We find the ground state of equation ([7j) by imaginary time evolution on a 3D grid [27J . 
The dipolar interaction part is efficiently computed in momentum space by means of fast 
Fourier transformations, as the mean field potential $dd has the form of a convolution. 

Although both the Bose-Hubbard Hamiltonian and the Gross-Pitaevskii equation 
are derived from the same multi-particle Hamiltonian in second quantized form there 
is no direct link between the two models. This implies that there can be no relation 
between a, add, N and the Bose-Hubbard parameters Uq, U\ since the latter depend on \1> 
which is not fixed in the mean-field approach. 

2.4- Stability 

Dipolar quantum gases have a complex stability behaviour [25] which leads to some 
peculiarities when treating them numerically. Typically, a critical scattering length a crit 
can be defined, which depends on the geometry of the external potential and the strength 
of the dipolar interaction [25]. For all scattering lengths a < a cr u the system is unstable 
and no ground state can be found. The crossing of the stability threshold leads to a 
collapse of the condensate wave function which can easily be identified in the simulation. 
The collapse in a single harmonic trap has been studied in detail |25j. Above the critical 
value of the scattering length, a stable solution can be found for any a > a cri f 
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As we are going to trace out the stability threshold of the triple-well configuration 
we need to assure that the simulation is not crossing any unstable regions during the 
imaginary time evolution. We proceed as follows: The simulation is set up with an 
initial Gaussian wave function which spreads over all wells. The spreading is such that 
the width in z-direction is equal to the spacing between the wells to assure a certain 
fraction of particles in each well. We stress that the final state does not depend on the 
chosen initial wave function. The Gaussian can even be placed asymmetrically over one 
of the outer wells and the imaginary time evolution still yields the same (symmetric) 
ground state. 

In the first sequence of imaginary time evolution all interactions are set to zero 
(a = a dd = 0) and we reach the ground state for the non-interacting case, see figure [3] 
for add = 0. This is our starting point to reach any point in the parameter space. The 
diagrams in figures H] and [5] are scanned line by line from right to left. We start at a high 
scattering length a with add still set to zero to assure that we are in the stable region. 
After the ground state for the contact interacting case is reached we select the final value 
of add and probe one horizontal line in the diagram by subsequent runs of imaginary 
time evolution while decreasing a until the wave function collapses to basically one grid 
point. At this point the simulation has reached the critical scattering length a crit . 

2.5. Geometry 

Performing the numerical simulations, we have to choose a fixed geometry and strength 
for the external potential, like the widths of a single well in all spacial directions 
w x ,w y , w z and the depth of the potential Vq. We reduce the parameter space by choosing 
reasonable values for the parameters, having symmetries as well as physical limitations 
in mind. 

The width in ^-direction is restricted, as the different wells are not clearly distinct 
for w z 3> 1/2 (remember that we are measuring lengths in units of the spacing between 
two wells). Oppositely, for w z <C 1/2 the tunneling is too low to reach the ground state 
of the system with imaginary time evolution (or in an experiment). In the simulations 
we will therefore set w z to a value of 1/2. We ask for two conditions to fix the values 
for the two remaining widths w x and w y . As we focus especially on inter-site effects, 
changing the polarization direction (from an attractive inter-site coupling to a repulsive 
one) should not change the on-site effects 0. To satisfy this requirement, the width 
in one of the directions perpendicular to z has to be equal to w z . Without loss of 
generality we define x to be the polarization direction for the "repulsive geometry" (z 
for the "attractive" case). Therefore we need to set w x = w z . 

The second condition concerns the stability. To see a large variety of states we 
want the stability of a single well to be higher than in the spherical case (lower critical 

| There might be changes in the on-site energy due to second order effects: if the changed inter-site 
coupling leads to a different shape of the on-site wave function. Strictly speaking, this condition can 
only be fulfilled for a single well. 
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scattering length a cr i t for the same add)- To fulfil this, the remaining width w y has to 
be larger than the two others [28]. This leads us to cigar-shaped traps which are placed 
side by side, as shown in figure [J |(b)| and (c) . We fix the trap aspect ratio to a value 



of oo z /ujy = 1/8 as this turns out to be a reasonable value for an experiment, too. We 
stress that simulations with different aspect ratios u z /u y < 1 do not show a qualitatively 
different behaviour. 

For the depth of the potential there are also certain limitations. If Vq is too low, 
the potential is not able to trap the particles. If it is too large, the tunneling rate is 
suppressed (see above). It turns out that V = 80 is a reasonable value which allows for a 
large diversity of ground states. Again, additional simulations show that the behaviour 
is not sensitive to the precise value of this parameter, even quantitatively. 



2.6. Interaction 

We simulate the dimensionless Gross-Pitaevskii equation (J7J). As all parameters of the 
external potential are fixed, there are only two free quantities. These are the values of 
the contact and dipolar interaction strength given by the dimensionless products aN and 
Odd-^- Note that it is not necessary to change the number of particles N independently. 
In the following we will present simulations where we change both aN and a^dN to scan 
the remaining parameter space. Note also that in our dimensionless units the values of 
a and Odd depend on the spacing between two lattice sites. 



3. Results 



To analyze the structure of the states we plot the ratio n = (n\ + n^)/N in analogy 
to [21]. In the simulation we calculate the occupation numbers n« by integrating the 
density n(r) over the volume of the i-th well. We have divided the whole volume of the 
simulation into three parts such that Yl n i = N. 

As we have a finite tunneling rate due to the finite potential depth and spacing, 




Figure 3. Purely dipolar states for a = and different values of addN in the repulsive 
case. Particles move to the outer wells as the dipolar inter-site repulsion grows. 
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the states found with the mean-field calculations are always symmetric (n% = n 3 ) with 
respect to the central well. The two asymmetric states in the D phase found in the 
Bose-Hubbard approach are only present for tunneling J = 0. For J > the symmetric 
and anti-symmetric combination of both states split in energy and yield a symmetric 
density distribution. 

Let us first discuss the non-interacting ground state of the triple-well system. In a 
simple 3-mode approach we can use localized wave functions = ^/njjNfa, centred at 
the z'-th well, as defined in section 12.11 If the ground state energy of a single well is E 
and the overlap integral for neighbouring wells is J, we have to diagonalise 

(e -J \ 

H= —J E —J , (9) 
V —J E J 

from which we immediately find the ground state (1/2, l/v2, 1/2) with occupation 
numbers n<i = N/2 and n\ = n% = N/4, giving a ratio of n = 1/2. Close to the origin 
of the diagrams in figure H] we find indeed states with n ~ 1/2 (see also figure [3] for 
a &&N = 0). Note that the non- interacting ground state in the Bose-Hubbard model for 
Uk = is given by 

(a\/V2 + 4/2 + a\/V2) N |000) 

which also yields (712) = N/2 and (ni) = (n 3 ) = N/4. 

Adding a repulsive contact interaction leads to a flattening of the density profile in 
the sense that we expect to have a uniform distribution n\ = n-i = 713 = N/3 (or a ratio 
of n = 2/3) for large scattering lengths. For a^d — and a > we observe states with 
1/2 < n < 2/3. Every state that is observed for add > which has a ratio n outside 
this interval is thus a clear indication of the dipolar inter-site effects. 

3.1. Repulsive inter-site interactions 

Figure H] shows an overview of the states found by imaginary time evolution for the 
geometry with repulsive inter-site interactions. We find the whole spectrum < n < 1. 
Once n approaches the value of or 1, the states become unstable. In particular, we 
observe ground states with a ratio of n > 2/3, implying that there are fewer particles 
in the middle well than in the outer ones (ri2 < n\ — n 3 ), a clear indication of the inter- 
site repulsion. States with n«0 appear even in the purely contact interacting case for 



G dd-^ = 0, see Figure 4(b) , and are therefore less suited to demonstrate the long-range 
nature of the interaction. The on-site attraction for negative a is enough to concentrate 
the atoms in the central well until the condensate finally collapses for aN —0.055. 



Figure 4(b) also reveals that a threshold value exists at add-^V ~ 0.3 with a sudden 
change of behaviour. For dipolar interactions weaker than this critical value, the ratio 
n is always lower than 2/3. This region is dominated by the on-site interactions. For 
decreasing contact interaction the ratio smoothly approaches n ~ and finally collapses. 
Contrarily, for an interaction strength larger than addN ~ 0.3, the ratio n increases 
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Figure 4. The quantity n = (ni 



inter-site interactions, see figure 1(b) 



■ fi3)/N is plotted for the geometry with repulsive 
(a) Black coloured areas indicate regions of 



instability where no ground state could be found. For the purely contact interacting 
case (add = 0, a > 0) we find an equally populated state with m ~ ~ n^. For 
increasing add the state transforms into the state with — 0, ri\ = 713 = N/2 (white). 
This state only appears close to the instability border, (b) Horizontal cuts through the 
phase diagram on the left are shown for different values of a^dN, as indicated by the 
labels. The behaviour changes qualitatively around the threshold value addN « 0.3. 
Note that only a narrow region of the phase diagram is shown in the cuts. 



when lowering the scattering length until it finally reaches n ~ 1 and collapses. The 
behaviour of the border, separating stable from unstable regions, is also different below 
and above the threshold. Below, the instability is triggered by a collapse in the middle 
well as it holds most of the particles. The system gets more stable for higher values of 
addN (critical scattering length decreases with growing dipolar interaction). Above the 
threshold, the instability is caused by the particle flow to the outer wells which finally 
leads to a collapse in either the left or the right well. This is a clear signature of the 
inter-site effects. In this regime, the system gets less stable for growing dipolar strength 
(a cr it increases for growing dipolar interaction). We remark that the threshold value is 
not depending on the aspect ratio of the single wells but changes with the depth of the 
potential. 



3.2. Attractive inter-site energy 



Figure 5(a) shows an overview of the states for the geometry with attractive inter-site 
interactions. We can immediately see that the inter-site attraction is destabilizing the 
system as it has a larger a crit for most values of add- We observe the whole spectrum 
of ratios n between the equally populated state with n = 2/3 and n — 0. Approaching 
the value of n — 0, the states become unstable. In this situation, the collapse is again 
initiated in the middle well. As we do not observe any states with n outside the range 
< n < 2/3, the attractive geometry is not suited to demonstrate the inter-site effects 
doubtlessly. 
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Figure 5. (a) The quantity n = (ni+nz) / N is plotted for the geometry with attractive 



inter-site interactions, see figure 1(c) Black colour again indicates unstable regions. 
Close to the instability, states with n « 0.5 form which resemble the non-interacting 
ground state, (b) Combined plot which shows the repulsive case (upper part) and the 
attractive case (lower part, flipped add axis). Both situations match exactly if add = 
but also for larger add the plot combines to a consistent picture. Differences in both 
parts are solely caused by inter-site effects. 



3.3. Combined picture 



Figure 5(b) shows a combined plot of both cases (repulsive and attractive inter-site 
interaction). For add = both situations are identical. Differences in the upper and 
lower part are solely caused by inter-site effects, as the geometry of the triple-well 
potential was designed in such a way (the on-site energy does not change when rotating 
the polarization direction). 

We now compare our findings to the results of the Bose-Hubbard approach. Within 
the mean-field theory and our numerical simulations we observe the counterparts of all 
states of phase A, as well as states close to those of the phases B and C. For phase 
D with its two degenerate states n\ = ri2 = N/2 and ri2 = = N/2, the comparison 
is a bit subtle. The mean ratio of both states is n = 1/2, which can be seen in the 
simulations. However, these states are related to the non-interacting state which also 
has n = 1/2. The mean- field approach is unable to distinguish both cases and finds 
only the symmetric states. 

As we do not observe extended regions in the phase diagram with n = or n = 1 
we conclude that the regions which correspond to the phases B, C and D are unstable, 
an aspect which is not observable in the Bose-Hubbard approach. We remark, however, 
that the presented theory is only valid for large particle numbers. For small samples 
with few particles, the Bose-Hubbard approach is well justified and stable phases should 
appear. 

Even in the absence of extended phases, the clear indication of inter-site effects 
is still visible. Especially the ground states with n > 2/3 above the threshold value 
at a^dN ~ 0.3 do not appear for purely contact interacting condensates and should 
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therefore be considered as a strong evidence for dipolar inter-site interactions. 

Finally, we justify the range of values used for the parameters aN and a^N. We 
recall that the unit of length we use is the spacing between two wells. We adopt the 
suggested experimental parameters of [21] where a spacing of I = 1.7|J.m is used. For 
2000 52 Cr atoms with a dipolar length of add ~ 0.79 nm, the (dimensionless) value of 
a ddN ~ 0.93 has the right order of magnitude. The scattering length a w 5.8 nm leads 
to the dimensionless value aN « 6.8. This value, however, can be tuned precisely by 
means of a Feshbach resonance [29], thus allowing for the detection of the interesting 
states close to the border of instability. Smaller samples of atoms could still provide 
the necessary dipolar interaction strength when using different species like Dy with a 
dipolar length of add ~ 7.1 nm [30J. 




4. Conclusions 

In this paper we studied the properties of a dipolar quantum gas, loaded into a triple-well 
potential. Using the numerical solution of the nonlocal Gross-Pitaevskii equation, the 
phase diagram of the system has been obtained. In particular, we identified the range 
of parameters where the ground states reveal strong inter-site effects and we traced 
out the instable regions in the phase diagram. We find that ultracold gases of atoms 
with a high magnetic moment, like Cr or Dy, are suited to demonstrate these features. 
The presented ground states were compared to the results of the dipolar Bose-Hubbard 
model. 
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